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Abstract 



The Escalator Boxcar Train (EBT) is a numerical method that is widely used in theoretical 
biology to investigate the dynamics of physiologically structured population models, i.e., models in 
which individuals differ by size or other physiological characteristics. The method was developed 
more than two decades ago, but has so far resisted attempts to give a formal proof of convergence. 
Using a modern framework of measure- valued solutions, we investigate the EBT method and show 
that the sequence of approximating solution measures generated by the EBT method converges 
weakly to the true solution measure under weak conditions on the growth rate, birth rate, and 
mortality rate. In rigorously establishing the convergence of the EBT method, our results pave 
the way for wider acceptance of the EBT method beyond theoretical biology and constitutes an 
. important step towards integration with established numerical schemes. 
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1 Introduction 

The population dynamics of ecological and biological systems are often described by an ordinary 
differential equation of the form 

where N = N(t) is the total population size at time t, /3(N) is the birth rate, and fi{N) is the 
mortality rate, both of which depends on the population size. The key assumption in this type of 
model is that every individual in the population is identical. This is clearly unreasonable in many 



situations, including cases where the gap between birth size and reproductive size is important. A 
more accurate description of the population dynamics can be given by physiologically structured 
population models (see e.g., [H]). In these models, the birth rates, death rates, and growth rates 
of individuals depend on their physiological state x G 0, where f2 is the set of admissible states. 
In general, these states can represent any aspects of individual physiology such as age, size, mass, 
height, or girth. For the purpose of this manuscript, we will work with a one-dimensional state space 
that we think of as representing individual size, but other interpretations are possible and, as we note 
in the concluding discussion, we expect that our results can easily be extended to higher-dimensional 
state manifolds. 
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In order to specify a physiologically structured population model, we need explicit representations 
for the mortality growth, and fecundity rates of individuals as well as the initial population structure. 
We assume that these rates are respectively on the form /j,(x,E t ), g(x,E t ), and j3{x,E t ), where x is 
the size (or more generally the state) of the individual and Et is the environment that individuals 
experiences at time t. The environment is a key factor in the formulation of physiologically structured 
population models and can, for example, represent the total amount of nutrient available at time t 
or the size-specific predation rate, see e.g., |181 [5]. While the environment is often low-dimensional, 
it could potentially be infinite-dimensional as would for example be the case for the shading profile 
in a forest. Finally, we assume that all new individuals have the same birth size x b . With these 
assumptions, one can show (see e.g., [5]) that the density u(x,t) of individuals of state x at time t is 
given by the first order, non-linear, non-local hyperbolic partial differential equations with non-local 
boundary condition 

d d 

— u(x, t) + Q^ {g( x , E t ) u(x, t)) = -fi(x, E t )u(x, t), (la) 
g(x b ,E t )u(x b ,t) = p(t,E t )u(t,t)(%, (lb) 



u(x,0) = u (x), (lc) 

in which we assume that x b < x < oo and t > 0. 

The first numerical method designed specifically for solving physiologically structured population 
models was the inventively named Escalator Boxcar Train (EBT) [4J. Rather than approximating 
the solution directly, it approximates the measure induced by the solution. Regardless of its uncon- 
ventional solution methodology, the EBT method is widely used by theoretical biologists (see e.g., 
[21 [131 032 [2D] ) • One of the reasons for the popularity of the EBT method can be ascribed to the 
simple biological interpretation of the components of the scheme: the state-space is partitioned into 
initial cohorts and, for the ith cohort, the EBT method tracks its size iVj(t) and the location of its 
centre of mass Xj(t) (see e.g., [5]). The solution measure d(t '■= u(t,x)dx is then approximated by 

TV 

dCtwdCf = Y, N ^ 5 Mt)> ( 2 ) 

i=B 

where 5 X is the Dirac measure concentrated at x. The dynamics of the functions iVj and Xi will 
be defined in Sect. [21 The boundary cohort corresponding to i = B is treated differentially from 
the other cohorts to account for newborn individuals. In the original formulation, this included 
terms correcting for changes in the average mass arising from the inflow of newborn individuals. For 
completeness, we consider the original definition of the boundary cohort in Sect. [H 

The convergence of the EBT method has remained an open question since the method was first 
introduced in 1988. The most successful analysis was performed by de Roos and Metz [6] in 1991. 
They studied how well the EBT method approximates integrals of the form j^ L ip{x)u{x,t) dx for 
smooth functions ip, assuming that cohorts are not internalized (see Sect. [2]). The result does 
not assert the convergence of the EBT method but rather, in the language used by de Roos and 
Metz, that the EBT method consistently approximates integrals of the solution to (TTJ). One reason 
for the lack of progress is that the usual analytical techniques for analyzing finite element and 
finite difference schemes are not immediately applicable to the measure-valued case. Over the last 
two decades, however, the theory of structured population models has been extensively developed 
[12\ O El El [TU1 [TT1 114j . and for the first time a full analysis of the EBT method is within our reach. 

The aim of this paper is to rigorously prove the convergence of the EBT method. We show that 
the EBT method converges under far weaker conditions on the growth, death and birth functions 
than the conditions assumed by de Roos and Metz [BJ. Our arguments build on recent theoretical 
developments by Gwiazda et al. [13] that extend the classical concept of weak solutions to measured- 
valued solutions. In the following section, we describe the EBT method in full detail, define weak 
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convergence of measures, and define weak solutions to the physiologically structured population model 
([1]). In Sect. [3] we prove the convergence of the EBT method with dynamics of the boundary cohort 
as introduced in this paper. Our convergence result is then extended to the original definition of 
the boundary cohort in Sect. 01 We conclude by placing our results into context and by highlighting 
promising directions for future work. Theorem 1141 and Theorem [16] are the main results of this paper. 



2 The Escalator Boxcar Train 



The EBT method is a numerical scheme for solving physiologically structured population models 
(PSPMs, see e.g., [IE]). While there are many possible formulations of PSPMs, several of which are 
described in the excellent book by Metz and Diekmann [18], we consider the numerical solution of the 
one-dimensional PSPM with a single birth state Xb defined by (jlap . (jlbj) . and (fTc]) , The EBT method 
determines an approximate measured-valued solution to the PSPM as a linear combination of 
Dirac measures, 

JV 
i=B 

Each of the terms in the approximation can be interpreted biologically as a cohort composed of N{ 
individuals with average individual state (e.g., size) Xi at time t. As individuals give rise to offspring 
with state Xb at birth, we need different definitions for internal cohorts and the boundary cohort. 

The internal cohorts are numbered i = B+l, N. These cohorts are chosen at time t = so that 
Cq^ converges weakly to the initial data uo(x)dx as N — > oo. This is always possible since finite linear 
combinations of Dirac measures are dense in the weak topology [1, Volume II, p. 214]. Thus, we need 
not restrict ourselves to initial data prescribed by a function uo(x), but can extend our analysis to 
general positive Radon measures vq. Without loss of generality, we will assume that the total mass 
Cq ([xi,, oo)) = zA)([xft, oo)) for all N. The boundary cohort is the cohort with the lowest index B. At 
time t = 0, B = and we assume that Nq(0) = and Xq(0) = Xb- As time progresses, additional 
cohorts with negative index will be created through the process of internalization described further 
below. 

The dynamics of the internal cohorts are given by 

^ = -v(X l ,( N )N l , (3a) 

f =9^,^, (3b) 

where we have assumed a direct dependence of the vital rates on the solution measure, ( N = ^ , to 
represent environmental feedback. Similarly, but in contrast to the original formulation of the EBT 
method by de Roos [3J, the dynamics of the boundary cohorts follow 

^ = -»(X B , ( N )N B + P(Xi, t N )N it (4a) 

i=B 

dXs , v j. _/y 



dt 



g(x B ,C), (4b) 



where the sum is taken over all cohorts including the boundary cohort. This sum reflects the offspring 
produced by the total population. In line with their biological interpretations, we henceforth assume 
that all vital rates, the mortality rate the fecundity rate f3, and the growth rate g, are non-negative. 

With the EBT method defined as above, both the width and the number of individuals in the 
boundary cohort will increase over time which eventually introduces an unacceptably large approx- 
imation error. For this reason, the boundary cohort must be internalized sufficiently often. This 
implies that the number of cohorts will increase following internalization. The new boundary cohort 
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is at the time t of the internalization given by Ng(t) = and Xs(t) = Xb, where B equals the 
index of the old boundary cohort decremented one step. At the same instant, the previous boundary 
cohort becomes an internal cohort. To prevent the number of internal cohorts from exceeding com- 
putationally acceptable bounds, internal cohorts may be removed when the number of individuals 
has declined sufficiently. Removal of internal cohorts is important for numerical implementation but 
will not be considered in this manuscript. 

The EBT method differs from traditional numerical schemes in that it aims to approximate the 
solution as a measure of point masses. Before we can discuss the convergence of the EBT method, 
it is necessary to extend the classical concept of a weak solution to measures. This extension builds 
on earlier work by Gwiazda et al. [2] (see also [31 [15]) and Chapter 8 of the monograph pQ. We 
will work with the cone all finite positive Radon measures denoted A4+(£l), where is a metric 
space consisting of all admissible individual states. In our presentation, we assume f2 = [xf,, oo] and 
we think of x £ 0, as the size of an individual. An important reason for working with finite Radon 
measures is that their behavior at infinity is tightly controlled: for each e > 0, there exists a compact 
set K e such that /i(Tl\K t ) < e. 

Since the EBT method approximates the true solution as a measure of point masses, the natural 
mode of convergence on A4 + (VL) is weak convergence^. 

Definition 1. A sequence of measures {nk} on O converges weakly to a measure if 



4>(x) dfi k (x) ->■ / 4>(x)d/i(x), 
a Jn 

as k — > oo for all bounded continuous real functions <fi on VL. 

The weak convergence defined above induces a topology associated with the Kantorovich- Rubinstein 
metric: 



p(p,v) = sup\ / 4>(x) d({i - v) 



n 



E C °°(R), U\\ w i,oo < 1 



in which ||^||^i,oo = ||0|lx,oo + H^'H^oo- This is also known as the flat metric. With this metric, 
A4 + (fl) is a complete metric space (see [21 Def. 2.5]). 

Analogously to weak convergence, we define weak continuity as follows: 

Definition 2. A mapping Q : R + — > A4 + (Q) is weakly continuous in time if for all bounded 
continuous real functions (f) on Vt, 

4>{x)dQ t , 



n 

is continuous in the classical sense as a function oft. 

With these two topological notions in place, we are in position to define measure- valued solutions 
to the PSPM ([I]): 

Definition 3. A mapping Q '■ [0,T] — > oo)) is a weak solution to ([1]) up to time T if Q is 

weakly continuous in time and 

OO POD 

(j)(x, T) d(x(jx) — / 4>(x,0) dvo(x) = 
b ' x b 

J (^^'^ + 9( x ^Ct)^(x,t) - fi(x,Ct)<l>(x,t)j d( t (x)dt 

pT poo 

+ / <f>{x h ,t) / /3(x',( t )d( t (x')d( t (x)dt, (5) 

J Jx b 

for all <f> E Co°(R + x [0,T]). Here, vq E M+{VL) is the initial data at time t = 0. 



1 There are two natural notions of convergence on M+(Q) — strong convergence and weak convergence. Strong 
convergence is unsuitable for our purposes as, for example, the sequence of Dirac measures 8\/ n does not converge to 
8q as n — > oo in the strong topology. 
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Remark 4. The definition above was inspired by Gwiazda et al. We differ in that we use 

smooth test functions, but note that these are dense in the space C 1 (~1 W 1 ' 00 used in \1J$ - 

Remark 5. The dependence on the environmental feedback variable E in (pQ) is represented here by 
a direct dependence on the solution measure Q. 

In order to show the convergence of the EBT method, we will recast the definition of a weak 
solution. Let < t\ < t 2 < T and v G M + (Q). For a given test function (f> € Cq°(R + x [0,T]) and a 
family of measures at, we define the residual 

poo poo 

R < f > (a t ,u,ti,t 2 ) = / <j>(x, t 2 ) da t2 (x) - / <f>(x, h) dv{x) (6) 



i 2 poo 



h Jx b 
t 2 



— (x,t) + g(x,( t ) — (x,t) - (j,(x,( t )(f>(x,t) ) da t (x)dt 



<p(x b ,t)^J /3(x',Ct)da t (x')j dt, 



where the measure v is interpreted as the initial data at time t = t±. Clearly, if R ( p(at,i / o,0,T) = 
for all test functions <f> and the family of measures at is weakly continuous in time, then at is a weak 
solution to ([T]). We will sometimes write RAat) meaning RAa t ,UQ,Q,T). 



3 Convergence of the Escalator Boxcar Train 

We establish the convergence of the EBT method in five steps: (1) At each fixed time t, the sequence 
of approximating EBT measures contains a subsequence which converges weakly to a positive Radon 
measure Ct- (2) We find a subsequence that for all t converges weakly to a mapping Q that is 
weakly continuous in time. (3) The residuals of the approximating EBT measures converges to 
the residual of Q for any test function. (4) The residual of the approximating EBT measures ^ 
converges to zero, and hence the measure Q is a weak solution. All that remains is then to show 
that the entire sequence of approximating EBT measures converges weakly to We do this by (5) 
assuming the existence of a unique weak solution to the structured population model and showing 
that a contradiction will otherwise result. In all the following lemmas, we assume that the birth 
rate, growth rate, and mortality rate are non-negative, bounded, and Lipschitz continuous functions 
of the individual size x. In addition, we need three assumption pertaining to the feedback from the 
population-level to individual vital rates: 

sup |/3 (x, a) - P(x,\)\ < C/3 p(a,X), 

X 

swp\g(x,a) - g(x,X)\ < C g p{a,X), 

X 

sup \n(x,a) - fj,(x, A) | < C M p(a, A). 

X 

The three requirements above assert Lipschitz continuity in A4+(0) equipped with the Kantorovich- 
Rubinstein metric. 

Lemma 6 (Step 1). For each t € [0, T], the sequence {Ct 1 } of approximating EBT measures contains a 
weakly convergent subsequence. In fact, any subsequence } of{^} contains a weakly convergent 
subsequence. 

Proof. By Prohorov's Theorem [TJ, it is enough to show that the sequence {C^} is uniformly bounded 
in the variation norm and is uniformly tight. As the measures are positive by construction, this 
amounts to showing that Ct ([xfc, oo)) is uniformly bounded in N, with limM->oo sup M (^((M, oo)) = 
0. An biological interpretation of these requirements, which we will build on in the proof, is that 
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the abundance and typical size of individuals in the population are bounded from above. Letting 
Pn(s) = (!? \[xb,oo)) it follows that 

N N N 

PM = £ Ni(s) = £ p(Xi, C?)Ni{s) - J> Cf) Ni{3) < 

i=B i=B i=B 

N N 

< Pi*, g)Ni(s) < /3 SU p W = ^u P Pn(s), 

i=B i=B 

where /3 sup is the supremum of /?, i.e., the maximum individual birth rate. The above inequality 
holds for all s G [0, T] except at the finite number of times, where boundary cohorts are internalized. 
At these points, the function P/v is continuous. Thus, < P/v(f) < Pjv(0) exp(/3 sup T). Hence 
Pjv(i) = Gt([%bi °°)) is uniformly bounded on [0,T], since Pjv(O) is independent of N. (Recall that 
in Sect. [3] we assumed that the initial mass should be independent of TV" and equal to that of the 
population measure given as initial condition.) 

To prove limM->-oo supjv (J^( (M, oo) ) = 0, we first show that the statement is true for t = 0. Let 
e > be given. Since the initial data vq is a positive Radon measure and thus tightly controlled at 
infinity, we may choose M\ large enough such that vq( (Mi, oo) ) < e/2. Pick any continuos function 
ip on [xb, oo) satisfying < (p(x) < 1 with tp(x) = 1 for x > M\ + 1 and (f{x) = for x < M\. Then 

/>00 /'OO 

Co ([Mi + 1, oo)) < / if dtg < Lpdv + e/2 < e, 

if we choose N > Nq for some sufficiently large Nq, since ^ converges weakly to uq as ./V — > oo. 
To account for the measures with N < iVo, we choose so large that ([M2, 00)) < e for N = 
1,2, ...,Nq. Finally, we choose M as the largest of the two numbers Mi + 1 and M2. 

To prove the statement for a general time t G [0, T], we first note that the center of mass 
and abundance at time t of any internal cohort i > with Xi(t) large enough can be estimated 
with their respective values at time t = 0. Specifically, -Xj(i) < -Xj(O) + tg sup , where g sup is the 
supremum of the growth rate g, and iVj(i) < A^(0). Combining these two estimates, we have that 
C/^(M, 00) < Co^(M — tg sup ,cxo) and the first assertion of the lemma follows. Finally we note that 
the above argument holds for any subsequence of {C/^}- This concludes the proof. □ 

Lemma 7 (Step 2). The approximating EBT sequence {Ct^} contains a subsequence which, for each 
t G [0, T], converges weakly to a positive finite measure £t- The mapping £ t : [0, T] — > A4 + (il) is 
weakly continuous in time. 

Proof. Let {qk}^ = i be an enumeration of the rational numbers in [0, T]. According to Lemma[6]there 
exists a convergent subsequence {C qi J } of {C^}- Repeating this argument, there exists a convergent 

N 2 N 1 N k 

subsequence {£ 92 J } of {Cg 2 J }• Proceeding by induction, we obtain for each k a sequence {Cj fc J } 
which converges weakly to Q qk and is a subsequence of all preceding sequences. Inspired by Cantor's 

diagonalization argument we define the sequence Q '■= Ct k ■ ^ follows that for each rational t € [0, T], 
this sequence converges weakly to a measure Ct- 

We will now show that the subsequence also converges to a positive finite Radon measure for 
all real t € [0,T]. We first show that for each fixed test function eft G C°°(M) and each time t, the 
sequence of real numbers 

poo 

/ <MCt fc , (7) 

converges as k — > 00. It then follows from classical results in the theory of distributions, e.g., |17[ 
Theorem 2.1.8 and Theorem 2.1.9], that Q converges weakly to a positive measure Q. This will turn 
out to be the desired measure. 
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To prove convergence of the sequence ([7]), we first note that for fixed k, the measure Q is weakly 
continuous in time since each iVj(.) and Xi(.) are continuous functions. Let t € [0, T] and <j> be a test 
function. Given e > we get 



4>dC t 



J-b 



< 



< 



peso roc roc roc roc roc 

/ <j> d(i - <f> dC* + / dQ - cP dC k q + / < fc - / 
iij Jxi, At jij «j 



for any j, k, and g. Noting that the birth rate and mortality rate are bounded, we can use the same 
argument as in the proof of Lemma [6] to show that the first and last term above are bounded by a 
constant multiple of \t — q\. In particular, this constant depends on neither j nor k. Choosing q as 
a rational number sufficiently close to t these two terms will be smaller than e/2. Finally, since q is 
rational, we may choose j and k large enough to make the middle term less than e/2. Thus, we have 
established the Cauchy property for the sequence ([7]), which hence converges for all test functions </>. 
This shows that Q converges weakly to a bounded positive Radon measure Ct for all t G [0,T]. 
Using the same idea as above, we see that Q is weakly continuous in time. Specifically, 



d( s 



*b 



*b 



< 



< 



dCs 



dC k s 



+ 



dC 



dll 



■i-b 



+ 



d& 



d(t 



X'f, 



where again the middle term is bounded by a constant multiple of \t — s\ independent of k. Finally, 
the first and last term can be made arbitrarily small as a consequence of the weak convergence of Q 
to ( s . □ 

Lemma 8. Assume that the sequence converges weakly to a finite Radon measure If ip G 
Cq°(M+ x [0, T]) then, for every bounded Lipschitz continuous function f satisfying 

sup|/(x,a) - f(x,X)\ < C f p(a, A), 

X 

for all a, A € A4 + (0,), we get 



T 



p(x,t)f(x,tf)dtf(x)dt 



p(x,t)f(x,(t)d( t (x)dt, 



as k tends to infinity. 
Proof. We have 



p(x,t)f(x,(t)d(Hx) 



j-b 



+ 



x b 
oo 



■I'b 



v(x,t)f(xXt)d( t k (x)+ 

p(x,t) (f(x,(t)-f(x,(t)) d$(x). 



(8) 



In the first term on the right hand side, the function p(x, t)f(x, Q) is bounded and Lipschitz contin- 
uous in x. Hence, it can be approximated by a sequence {p m } of functions in C^°(M + x [0,T]) that 
converges pointwise and in W 1,co -noi:m. Such a sequence can, for example, be constructed through 
convolution. As the first term would vanish if \\<p(-, t)f(-, Q)\\yyi,<x, = 0, we can assume that this is 
not the case. We then get 



<p(x,t)f(x,£t)d$(x) 



xb 



poo 

/ lim p m (x,t) dCf(x) -- 

Jx h m ^°° 

ip m (x,t)d(t(x), 



lim 

m— >oo 



■i'b 
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where we have used Lebesgue's dominated convergence theorem. Hence, 



<p(x,t)f(x,C t )d${x) 



■l-b 



<p(x,t)f(x,(t)d( t (x) 



lim 

m— >oo 



V m {x,t)d{$ -&){*) 



x b 



J-'b 



< 



< lim \\ip m (., t)\\ w i >00 p(C t ~ Ct) = 
m— >oo 

= \H;t)f(.Xt)\\ W ^p(tf-Ct)^0, 

as k tends to infinity. Thus the first term converges to 

(p(x,t)f(x,(t)d(t(x). 

x b 

It remains to show that the second term in §S§ vanishes as k — > oo, 

ip(x,t) (/(x,C t fc )-/(x,Ct)) d$(x) 
<p{x,t) (f(xXt)-f(xXt)) |Ct fc (K,oo)) < 
< sup|y>(x,t)|sup f(x, Ct) ~ f(x, Ct) Ct fc (N,oo)) < 

X X 

<C v C fP (Ct Ct)CH[*b,oo)). 

Since Ct converges weakly to Ct, it follows from Gwiazda et al. [21 Theorem 2.7] that Ct([ x bi°°)) i s 
uniformly bounded and p(Ct,Ct) tends to zero as k tends to infinity. Since the above calculation is 
done pointwise in t, the lemma follows from Lebesgue's dominated convergence theorem. □ 

Lemma 9 (Step 3). Assume that the sequence Ct converges weakly to the finite Radon measure Ct- 
Then the residual RACt) converges to RACt) f or a ^ test functions (j) £ Co°(^+ x [O,? 1 ])- 

Proof. Consider 



< sup 

x 



roo poo 

R<p{Ct)= <Kx,T)d${x)- / <P(x,0)dv Q (x) 

Jx b Jx b 



(9) 



T i-oo 
Jx b 



-l(x,t)+g(x,Ch £(x,t)-p(x,tf)Ht,x)) dtf(x)dt 



+ 



£ (j>(x b , t) i3(x', Cl ) d&V)) dt = I-II-III + IV. 



The first term converges by definition of weak convergence and the second term is unchanged. The 
third and fourth term converge by Lemma [H □ 

Lemma 10. Let < t\ < ti < T and v E Ai+(£l). Assuming that no internalization is done in the 
interval (ii,^), then for any test function <fi we have that 



N 



+ 



(x, t\) dv(x)+ 



j-b 



N 



(<t>(x B (t),t) -<Kx b ,t))J2p(Xi(t),c?) Ni(t)dt, 

i=B 



where the sum is taken over all cohorts, including the boundary cohort. 
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Proof. We write the residual ^ as 

roo roo 

R 4 ((?,v,t 1 ,t 2 ) = <t>(x,t 2 )d(g(x)- <f>(x,h)du(x 

JXh JXh 



x b Jx b 

1 / {M^t)+9(x,C t N ) 0i(M)-/*(*>Cf)0OM)) d( t N (x)dt 

ti J x b 
t-2 / roo 



t\ \Jx b J 

= I(( t N 2 )-II(v)-III(( t N )-IV(C t N )- 

Here we have used the shorthand notation 4>i(£,t) = d(j)(£,t)/dx and fai&t) = d(j>(£,t)/dt. 
Recalling that 

N 

Cf = Z>i(^(t)> 

i=B 

we get 

N 



j(O = E^^)0(^(t 2 ),t 2 ), 



i=B 



ni(C t N ) = X: / mt) (<t>2(Xi(t),t) + g(Xi(t),(t N ) MM*),*) - i-i(xi(t),Ct N )HXi(t),t)) dt 

i=B Jtl 

N 

= ni B (( t N )+ HI ^)- 

i=B+l 

Now, by (|3J), we have 

= j t [Niit^X^^dt = Ni(t2)<KXi(t2),t 2 ) - JV<(ti)0(Xi(ti),ti). 

Thus 

N N 

HO- E u h{C?) = N B {t 2 )^{X B {t 2 )M)+ E NiihWXiit&h). 

i=B+l i=B+l 

In the same way, but now also using Q, we get 

iii B ttt N ) = / - (JV B (t) <Kx B (t),t)) - <KXB(t),t)^p(Xi(t),c?) Ni(t)dt = 

J h at i=B 

= N B (t 2 ) HX B (t 2 ),t 2 ) - N B (h) 4>{XB{ti)M)- 

*2 



tl 



<A(X B (t),t)E/5(^(i),Cf) Ni(t)dt. 



i=B 



Since 



JV(Cf) = / <j>(x b ,t)J2P(Xi(t),(?) Ni(t)dt, 

Jtl i=B 



CONVERGENCE OF THE ESCALATOR BOXCAR TRAIN 



10 



we have 

-III B ((?)-IV(( t N ) = -N B (t 2 ) <j>{X B {t 2 )M +N B (tx) <f>(X b (h),tx)+ 



,t 2 n 
+ / {<j ) {X B {t),t)-4>{x b ,t))Y J P{X i {t),^) Ni(t)dt. 

Jtl i=B 



Summing up the calculations above, we get 

N 



i(0 - E in ^ N ) - nHcn - ma 

i=B+l 

N 

= N B (t 2 )<t>(X B (t 2 ),t 2 )+ Ni(h)HXi(-ti),tx)-N B (t 2 ) <j>{X B {t 2 ),t 2 )+ 

i=B+l 

rt 2 N 

+N B (h) </>(X B (tx),tx)+ / (<f>{X B {t),t) - <f>(x b ,t))J2KXi(t)Xn Ni{t)dt = 

Jtl i=B 
N t2 N 

Y J N i {ti)4>{X i {h)M)+ / (4>(X B (t),t) - <p(x b ,t))Y,P(Xi(t),(?) Ni(t)dt. 



i=B " 11 i=B 

Finally we get that 

v 



A " poo 

R^iC? ,u,tx,t 2 ) = ^N l (t 1 )0(X i (t 1 ),t 1 ) - / <j)(x,tx)du(x)+ 

i=B Jx b 

fh N 

+ / (<i>(x B (t),t) - d>(x b ,t))Y,p(Xi(t),c?) m)dt. 

Jtl ,= R 



□ 



Remark 11. The residual can be interpreted as the sum of the error arising from the discretization of 
the initial data and the error arising from the boundary cohort. In the interior of the individual state 
space, the EBT method gives an exact solution, i.e., there are no errors arising from the transportation 
of the interior cohorts. 

Lemma 12 (Step 4). With defined by the EBT method with internalizations at times ti = iT/n, 
we have that 

i^(C t Vo,0,T)^0, 
as N and n tends to infinity. Here v$ is the initial data at time t = to = 0. 

Proof. We first write 

n-1 

^(Ct^OjO,! 1 ) = R^C^ ,u Q ,0,ti) + Rcj, (C^ , Ct! ,U,tj+x)- 

1=1 

By Lemma [10] we have, 

^(Ci Y ,^o,0,ti) = ^iV J (0) ( / ) (x i (0),0)- / <j>(x,0)dvo(x)+ 

i=B 



+ / (^(X B (t),t) - c/ ) (x b ,t))Y J (3(Xi(t),^) Ni{t)dt, 



i=B 
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and 



A straightforward estimate now gives 

N 

|^(Cf,^o,0,T)| < ' 



ly poo 

^ ^(0)0(^(0), 0) - / ct>(x,Q)du (O + 

i=B Jx b 

+ Y, I \^ x B(t),t) -<Kx b ,t)\ ^P( Xj (t),( t N ) Nj(t)dt. 



i=0 H j=B 

The first term tends to zero by assumption as the number of initial cohorts, N, tends to infinity. 
Noting that x b = Xs(ti) and using that the growth rate is bounded, we get 



\<j>(X B (t),t)-<j>{x b ,t)\ < CtlXsit) - x b \ < Cte\t-U 



Hence, 



re— 1 P t- , , N 
i=0 tl j=B 



|0(X B (i),i)-^ 6 ,t)|^/3(^(t),Cf) Nj(t)dt < 

j=B 

n-1 

^ ] C<f>g \ti+l ~ ti\ C/3 UQ , 



i=0 

for the constant (Tg^ = /3 sup fo([xb, oo)) exp(/? sup T). Thus, the last sum is bounded by C(T)/n which 
also tends to zero as the number of internalizations tends to infinity. □ 

Remark 13. Examining the proof above, we see that the residual tends to zero whenever the maximal 
time between two internalizations of the boundary cohort tends to zero. Hence, we can relax the 
assumption that the times at which the boundary cohort is internalized are evenly distributed. 

Recalling that the initial cohorts are chosen to converge weakly to the initial data, we are now 
able to prove convergence of the Escalator Boxcar Train: 

Theorem 14. Assume that the assumptions on the birth, growth, and mortality rates in the beginning 
of Sect. \^hold. If the structured population model given by (jlap . (|lb|) . and (fTcj) has a unique solution 
(t, then the the solutions given by the EBT method converge weakly to Qt «s the number of initial 
cohorts tends to infinity and the maximal time between two boundary cohort internalizations tends to 
zero. 

Proof. (Step 5) We assume that the entire sequence £ t does not converge to Q. Then, in the weak 
topology, there exists an open neighborhood U of ( t , and a subsequence (f^ k of such that £ U 
for all JVjfc. From Lemma 101 1121 we conclude that {Ct fe } contains a convergent sub-sequence with a 
limit point not equal to Q, which is a contradiction since it would imply that the solution to the 
PSPM is not unique. □ 

The proof of convergence assumed exact solutions to the ordinary differential equations (ODEs) 
underlying the EBT method. In practical implementations, these need to be solved numerically 
which introduces small but finite approximation errors. We now extend the convergence proof to 
account for errors introduced by the underlying ODE solver. 

The following lemma is an immediate consequence of Lemma 

Lemma 15. Assume that C t ' = £i=B JV i l (*)< y xfc(t)- Ufor each t we have that Nf(t) -> N^t) and 
X?(t) -»• Xi(t) ash\0 then < h , u ,0, T) -»• i^(Cf, vq, 0, T) ash\ 0. 
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Combining the lemma above with Theorem [TJ] we finally have 

Theorem 16. Assume that the assumptions on the birth, growth, and mortality rates in the beginning 
of Sect. \3[hold. If the structured population model given by (jlap . (jlbj) . and (fTc|) has a unique solution 
Ct, then the the solutions C ' , given by the numerical integration of the EBT method, converges weakly 
to (t if the number of initial cohorts tends to infinity and the maximal time between two boundary 
cohort internalizations tends to zero, while h tends to zero sufficiently fast. 



4 The original definition of the boundary cohort 

Our study of convergence of the Escalator Boxcar Train in Sect. [3] assumed different dynamics of the 
boundary cohorts than was used in the original formulation of the method by de Roos We based 
our work on the assumption that the boundary cohort differed from the interior cohorts only in the 
addition of a term for the inflow of newborns. In this section, we consider the convergence of the 
EBT method under the original definition of the boundary cohort dynamics. 

While we simply assumed a dynamical system for the boundary cohort, de Roos formally derived 
the underlying equations. Consequently, the original dynamics for the boundary cohort reflect the 
reduction in center of mass that in reality accompanies an inflow of newborns. Moreover, as the 
center of mass is not defined as a physical quantity for an empty cohort, the equations were derived 
through series expansion around the size at birth. Thus, rather than tracking the center of mass 
X B (t) directly, de Roos considered a quantity tt b which roughly represents the cumulative amount 
by which the individuals in the boundary cohort exceed their birth size. This quantity is mapped 
onto the center of mass through the non-linear transformation 

( 7f B 

X B = \ N~ B +X - if7rB> °< (10) 
[ Xb, otherwise. 

The specific equations used for defining the boundary cohort were 

— — = -/i(x ft ,C )N B -K- tt b + 2^f3{Xi,C )Ni, (11) 

X i=B 

ditB i rN . AT . dg{x b ,( N ) N 

~Jf = 9K x b,C )N B ^ ir B -n(x b ,( )ttb, (12) 

with initial conditions N B = ir B = 0. We will assume that these are non-negative, as this is a 
natural requirement which can easily be enforced by an ODE solver if necessary. The appearance of 
partial derivatives in the expressions above, arising from series expansion around the size at birth, in 
conjunction with the non-linear transformation mapping tt b and N B onto X B , pose new challenges 
for proving convergence. As we will show, however, our proof of convergence can be tailored to 
accompany also the original definition of the boundary cohort. 

Note first that the only parts in the proof of convergence in which the equations defining the 
boundary cohort are used is Lemma [9] and implicitly in Theorem [T6j It therefore suffices to give new 
proofs of these statements. To this end, we require an additional lemma concerning the behavior of 
the quotient ir B /N B : 

Lemma 17. With N B and n B defined by (|lip and (|12p . we get 

0<X B -x b < Ct, 



for t € [io>^o + h] and some positive constants C and h which only depend on g, dg/dx and dfi/dx. 
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Proof. Prom the definitions of Nb and 7Tb we have 



— X— — — 

dt b ~ alt N B N B dt 



1 dir B ttb dN B 



N B dt 



9 + 



og kb 
dxN B 



7TB 7TB d{l 7r| 

^N B + ^N B + dxNl 



7T B 



N 



N 2 E« 



~ 9+ dxN B + dxN B n\ 2 ^ 
~ 9 dxN B dx \N B 

Remembering that X B = tt b /N b + a;^ we thus have X' B < a + 6(Xb — xj) + c(Jg — x b ) 2 for 
some positive constants a, b and c. Hence X' B < 2a when X B < X B for some positive X B . Since 
X B (0) = x b it follows that X B (t) < x b + 2at for i G [0, X|,/2a]. □ 

We now use this to show that 

Lemma 18. Assume that a new boundary cohort is created at time t = t±. For t 2 > t\ sufficiently 
close to t\, we have for all t G [ii , ^2] that 



N 



< 



N B (t) 



dX B (t) N , 
-g(X B (t),Q : 



dt 



and 



<c 2 (t 2 -t 1 ). 



(13) 



(14) 



Proof. Since N B is bounded, ir B = N B (X B — x b ), it follows from the above proof that |7rs(t)| < 
C{t 2 — t\) for some positive constant C. Hence, since also d/j,(X B (t)Xt*)/® x 1S bounded by the 
assumptions in [3], the statement (|14[) follows trivially. To show the first part of the assertion, we 
note that 



\ dg dfi 7T B 7T B A 

) = d-x 7 < B+ d-x 7TB N B -N- B ^ m - 

' i=B 



Since ir B and ir B /N B = X B — x b both increases at most linearly from zero, the assertion (fT3j) 
follows. □ 

The two lemmas above will be used to bound the residual between two internalizations. 

Lemma 19. Let < t\ < t 2 < T and v € M+(£l). For a given test function <j) G Co°(^+ x [0,T]) 
and a family of measures at- Assuming that no internalization is done in the interval (ti,t 2 ), then 

R<t>{Ct ',",h,t 3 ) = y ^2N i {t 1 ) ( j)(X i (t 1 ),t 1 ) - / </>(x,t 1 )dv(x)+ 



+ 



N 



(<Kx B (t),t) - </>(x b , t)) ^2 KM*), C t N ) W)dt+ 

dX B (t) 



i=B 



dt 



g(x B (t)X, 



N B (t) 

+ {^(X B (t),(?)K B (t)) 4>{X B {t),t)dt, 
where the sums are taken over all cohorts, including the boundary cohort. 



Mx B (t),t)+ 
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Proof. Examining the proof of Lemma [TU] we see that the boundary cohort only appears in the term 
Mb, 

Mb{C?)= f 2 N B (t) (MX B (t),t)+g(X B (t),Ct N ) MX B (t),t)- 
Jti 

-f,(x B (t),( t N )HXB(t),t))dt. 

This term is shown to be equivalent with 

*2 A N 

- (N B (t) <f>(X B (t),t)) - <j>(X B (t),t)J2HXi(t),Ct N ) Ni(t)dt. 

41 i=B 

Using the original definition for the boundary cohort dynamics, ([lip and (|12[) . we derive the required 
correction term 

t 2 a N 

j(N B (t) <t>(x B (t),t))-<t>(X B (t),t)J2HMt),( t N ) Ni(t)dt-III B (C?) = 

1 i=B 

t2 N B (t) (^^±-g( XB (t),(t N ) ) MX B (t),t) + 

+ f x 1 (X B (t),( t N )ir B (t)<f>(X B (t),t)dt. 

□ 

By Lemma \TEl we see that the correction term above is bounded by C(t2 — ti) 2 . Analogous to 
Lemma [T2[ we then have 

Lemma 20. With defined by the EBT method with internalizations at times ti = iT/n , we have 
that 

i^(CiVo,0,T) ->0, 
as N and n tends to infinity. Here is the initial data at time t = to = 0. 

The original definition of the boundary cohorts might prove more challenging from a numerical 
perspective. However, if we can determine numerically solutions to the equations of the EBT 
method such that the center of mass, Xg, now determined by the non linear transformation (jlOp 
converges to its true value, X B , as the step length h \ 0, the residual still tends to zero according 
to Lemma [T5l Hence, the numerical convergence follows as before. 



5 Discussion 



Enhanced biological realism and predictive ability of theoretical investigations are gaining importance 
as anthropogenic impacts are fundamentally altering the native environment of many organisms. 
Physiologically structured population models (PSPMs) are increasingly used to model and analyze 
biological systems. As these models account for the physiological development of individuals, they 
are better able to predict system dynamics. In contrast to simple unstructured population models 
such as the classical Lotka-Volterra equations, PSPMs often defy analytical investigations due to the 
non-local dependencies. There is thus a mounting need for numerical methods that can effectively 
uncover the underlying dynamics. The Escalator Boxcar Train (EBT) has been specifically designed 
for PSPMs and has three major advantages: it prevents numerical diffusion, it is relatively easy to 
implement, and the underlying equations allow for a natural biological interpretation. The method 
was developed more than two decades ago and has been used to study PSPMs ever since, but the 
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fact that convergence has never been formally proved might well have hampered its wider acceptance 
beyond the domains of theoretical biology. 

In this paper we have given the first rigorous proof of convergence for the EBT method. Our 
proof is given in a modern setting of measure- valued solutions (see e.g., |14|). This contrasts with 
previous efforts by de Roos and Metz [6] that were carried out in a classical setting and thus required 
additional smoothness assumptions. While their efforts fell short of proving the full convergence of 
the EBT method, the authors succeeded in showing that the method consistently approximates the 
true solution, i.e., that the local approximation error as measured through an arbitrary (but smooth) 
functional of the solution is bounded and vanishes in the limit of infinitely fine discretization of the 
individual state space. 

There are many possible extensions of the work presented here. A straightforward extension is 
to write down the corresponding proof for a higher-dimensional state space but with a single birth 
state. We believe that with more tedious calculations, one could prove the convergence also for the 
case of stochastic birth state. A more challenging extension is to consider stochasticity in individual 
development. On the population- level, this roughly amounts to diffusion and it is difficult to see how 
the EBT method should best be adapted to deal with this situation. Here, some inspiration might 
come from moving-mesh discontinuous Galerkin methods which, at least at first glance, appear to 
have similarities with the EBT method. A further extension is to consider different formulations of 
the boundary cohort. We initially proved convergence when the boundary cohort differed only by the 
addition of a fecundity term. While this works mathematically, it is natural to account for the fact 
that newborn individuals reduces the average size of individuals in the boundary cohort. The original 
formulation of the EBT method does account for this through a different definition of the boundary 
cohort, and as a second step we analyzed and proved convergence for this case. We believe that our 
proof can be extended to show convergence also for other formulations of the boundary cohort, as 
long as the flux of individuals is preserved. Analyzing convergence rates for different definitions of 
boundary cohorts would be an interesting extension of the work presented here. In particular, we 
believe that the series expansion around the size at birth underlying the original derivation of the 
boundary cohorts is not required, and that a direct evaluation at the center of mass might lead to 
even faster convergence. This could well be part of a more broadly encompassing study that explores 
convergence rates under different smoothness assumptions. A final important extension would be 
to consider vital rates that depend on the entire history of the population state up to the current 
time, rather than merely the current population state, as this would encompass cases with dynamic 
environmental feedback variables. 

Given the long tradition of partial differential equations (PDEs) in the physical sciences, it is not 
surprising that PSPMs were initially studied using this formalism. Efforts in the last decades have 
revealed, however, that the PDE formalism is not well-suited for considering questions of existence, 
uniqueness, and stability. For this reason, the cumulative formulation of structured population models 
[lU |9] was developed. It had the drawback, however, that a principle of linearized stability and the 
Hopf bifurcation theorem proved hard to establish [16]. Currently, it appears that renewal equations 
are well-suited for studying PSPMs [TOj, El [Til 05] • The work presented here has been developed 
from the PDE setting. We believe, however, that renewal equations are a promising framework 
for developing and analyzing numerical methods for PSPMs. A first step would be to recast the 
EBT method in this setting, after which the extensions outlined above could be considered. With 
interest in PSPMs now mounting, a historical opportunity exists for bridging biological theory and 
computational mathematics through the development of modern numerical methods for the 21st 
century. 
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